# 代码来源: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9470838/
#### 细胞类型自动化注释 ####
# 数据来源：http://biocc.hrbmu.edu.cn/CellMarker/download.jsp

library(dplyr)
library(openxlsx)

#### 1. 加载marker数据库 ####
load("./Human_cell_markers.RData")
cellmarker.file <- "./multi/multi_marker_genes_tsne_20PC.txt"
sampleid <- "multi"
# 读入单细胞分析中输出的cell marker 基因文件
marker.gene <- read.table(cellmarker.file, header = T, stringsAsFactors = F, sep = "\t")
marker.gene.sig <- marker.gene %>% filter(as.numeric(p_val_adj) <= 0.05)

cat("Totally marker genes:", length(unique(marker.gene.sig$gene)))
table(marker.gene.sig$cluster)

#### 2. 定位数据库中存在的基因 ####
marker.gene.sig %>% filter(gene %in% cell.markers.tb$geneSymbol) -> marker.gene.sel
marker.gene.sel$cellMarker <- apply(marker.gene.sel, 1, function(x) paste(unique(cell.markers.tb[cell.markers.tb$geneSymbol == x["gene"], 1]), collapse = ","))
cat("Totally find", length(unique(marker.gene.sel$gene)), "/", length(unique(marker.gene.sig$gene)), "genes in cellMarker db")

#### 3. 将marker基因与细胞类型结果写出到文件 ####
write.table(marker.gene.sel, file = paste0("DEG_marker_cells_tsne.txt"), row.names = T, col.names = T, sep = "\t", quote = F)

################################### 04.SingleR R包注释细胞类型###################################
counts <- pbmc@assays$RNA@counts
clusters <- pbmc@meta.data$seurat_clusters
ann <- pbmc@meta.data$orig.ident
# ref=get(load("ref_Human_all.RData"))
ref <- celldex::HumanPrimaryCellAtlasData()
singler <- SingleR(
    test = counts, ref = ref,
    labels = ref$label.main, clusters = clusters
)
clusterAnn <- as.data.frame(singler)
clusterAnn <- cbind(id = row.names(clusterAnn), clusterAnn)
clusterAnn <- clusterAnn[, c("id", "labels")]
write.table(clusterAnn, file = "singleR.clusterAnn.txt", quote = F, sep = "\t", row.names = F)
singler2 <- SingleR(
    test = counts, ref = ref,
    labels = ref$label.main
)
cellAnn <- as.data.frame(singler2)
cellAnn <- cbind(id = row.names(cellAnn), cellAnn)
cellAnn <- cellAnn[, c("id", "labels")]
write.table(cellAnn, file = "singleR.cellAnn.txt", quote = F, sep = "\t", row.names = F)



###################################05.monocle R包细胞轨迹分析###################################
#准备细胞轨迹分析需要的文件
monocle.matrix=as.matrix(pbmc2@assays$RNA@data)
monocle.sample=pbmc2@meta.data
monocle.geneAnn=data.frame(gene_short_name = row.names(monocle.matrix), row.names = row.names(monocle.matrix))
monocle.clusterAnn=clusterAnn
monocle.markers=sig.markers

#将Seurat结果转换为monocle需要的细胞矩阵，细胞注释表格和基因注释表格
data <- as(as.matrix(monocle.matrix), 'sparseMatrix')
pd<-new("AnnotatedDataFrame", data = monocle.sample)
fd<-new("AnnotatedDataFrame", data = monocle.geneAnn)
cds <- newCellDataSet(data, phenoData = pd, featureData = fd)
names(pData(cds))[names(pData(cds))=="seurat_clusters"]="Cluster"
pData(cds)[,"Cluster"]=paste0("cluster",pData(cds)[,"Cluster"])

#添加细胞聚类数据
clusterAnn=as.character(monocle.clusterAnn[,2])
names(clusterAnn)=paste0("cluster",monocle.clusterAnn[,1])
pData(cds)$cell_type2 <- plyr::revalue(as.character(pData(cds)$Cluster),clusterAnn)

#细胞轨迹分析流程
cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)
cds <- setOrderingFilter(cds, as.vector(sig.markers$gene))
#plot_ordering_genes(cds)
cds <- reduceDimension(cds, max_components = 2, reduction_method = 'DDRTree')
cds <- orderCells(cds)
#保存树枝的细胞轨迹图
pdf(file="05.trajectory.State.pdf",width=6.5,height=6)
plot_cell_trajectory(cds,color_by = "State")
dev.off()
#保存时间的细胞轨迹图
pdf(file="05.trajectory.Pseudotime.pdf",width=6.5,height=6)
plot_cell_trajectory(cds,color_by = "Pseudotime")
dev.off()
#保存细胞名称的细胞轨迹图
pdf(file="05.trajectory.cellType.pdf",width=6.5,height=6)
plot_cell_trajectory(cds,color_by = "cell_type2")
dev.off()
#保存聚类的细胞轨迹图
pdf(file="05.trajectory.cluster.pdf",width=6.5,height=6)
plot_cell_trajectory(cds, color_by = "Cluster")
dev.off()

#细胞轨迹差异分析
groups=subset(pData(cds),select='State')
pbmc=AddMetaData(object=pbmc, metadata=groups, col.name="group")
geneList=list()
for(i in levels(factor(groups$State))){
	pbmc.markers=FindMarkers(pbmc, ident.1 = i, group.by = 'group')
	sig.markers=pbmc.markers[(abs(as.numeric(as.vector(pbmc.markers$avg_log2FC)))>logFCfilter & as.numeric(as.vector(pbmc.markers$p_val_adj))<adjPvalFilter),]
	sig.markers=cbind(Gene=row.names(sig.markers), sig.markers)
	write.table(sig.markers,file=paste0("05.monocleDiff.", i, ".txt"),sep="\t",row.names=F,quote=F)
	geneList[[i]]=row.names(sig.markers)
}
#保存交集基因
unionGenes=Reduce(union,geneList)
write.table(file="05.monocleDiff.union.txt",unionGenes,sep="\t",quote=F,col.names=F,row.names=F)




pbmc <- RenameIdents(
  object = pbmc,
  "0" = "CD4+ cytotoxic T cell",
  "1" = "Kupffer cell",
  "2" = "Liver bud hepatic cell",
  "3" = "CD8+ T cell",
  "4" = "Endothelial cell",
  "5" = "Monocyte",
  "6" = "Liver bud hepatic cell",
  "7" = "CD8+ T cell",
  "8" = "Myofibroblast",
  "9" = "CD4+ cytotoxic T cell",
  "10" = "Hepatocyte",
  "11" = "Exhausted CD4+ T cell",
  "12" = "B cell",
  "13" = "Hepatocyte",
  "14" = "B cell",
  "15" = "Dendritic cell",
  "16" = "Hepatocyte",
  "17" = "Exhausted CD8+ T cell",
  "18" = "Cancer stem cell",
  "19" = "Liver bud hepatic cell",
  "20" = "Endothelial cell",
  "21" = "Liver bud hepatic cell",
  "22" = "Hepatocyte",
  "23" = "Kupffer cell",
  "24" = "Liver bud hepatic cell"
)


pbmc$celltype.group <- paste(Idents(pbmc), pbmc$TimePoints2, sep = "_")

pdf(paste0("CellCluster-slipt-time-tsne","PC.pdf"),width = 12,height = 15)
DimPlot(pbmc, 
        split.by ="TimePoints2", 
       reduction = "tsne",label = F,raster=FALSE)
dev.off()

pdf(paste0("CellCluster-slipt-time-umap","PC.pdf"),width = 12,height = 15)
DimPlot(pbmc, 
        split.by ="TimePoints2", 
       reduction = "umap",label = F,raster=FALSE)
dev.off()

pdf(paste0("CellClusterAll-group.pdf"),width = 20,height =8)
DimPlot(object = pbmc, 
        group.by="TimePoints2", 
        pt.size=0.5,reduction = "tsne",raster=FALSE)
dev.off()

pdf(paste0("CellClusterAll-group-umap.pdf"),width = 12,height =8)
DimPlot(object = pbmc, 
        group.by="TimePoints2", 
        pt.size=0.5,reduction = "umap",raster=FALSE)
dev.off()

pdf(paste0("CellClusterAll-group-tsne.pdf"),width = 12,height =8)
DimPlot(object = pbmc, 
        group.by="TimePoints2", 
        pt.size=0.5,reduction = "tsne",raster=FALSE)
dev.off()

pdf(paste0("CellCluster-sliptumap3.pdf"),width = 16,height = 8)
DimPlot(pbmc, 
        split.by ="TimePoints2", 
       reduction = "umap",label = F,raster=FALSE)
dev.off()

pdf(paste0("CellCluster-slipttsne3.pdf"),width = 16,height = 8)
DimPlot(pbmc, 
        split.by ="TimePoints2", 
       reduction = "tsne",label = F,raster=FALSE)
dev.off()

pdf(paste0("CellClusterAll-group4.pdf"),width = 12,height =10)
DimPlot(object = pbmc, 
        group.by="TimePoints2", 
        pt.size=0.5,reduction = "tsne",label=pbmc$celltype,raster=FALSE)
dev.off()

FeaturePlot(pbmc,feature ="INCENP",reduction = "umap",label = F,raster=FALSE)
VlnPlot(pbmc,feature = "INCENP",pt.size = 0)

INCENP
NEK11
AURKB
CCNA2
PRKAA2
CHEK1
CDK5
PRKCD
CDK1

vargene<-pbmc[["RNA"]]@var.features
write.table(file="vargene.txt",vargene,sep="\t",quote=F,col.names=F,row.names=F)

RNAgene<-pbmc@assays$RNA


#将seurat对象转换成monocle对象
##############################
data <- pbmc2[['RNA']]@data
pd <- new('AnnotatedDataFrame', data = pbmc2@meta.data)
fData <- data.frame(gene_short_name = row.names(data), row.names = row.names(data))
fd <- new('AnnotatedDataFrame', data = fData)
HSMM <- newCellDataSet(data,
                       phenoData = pd,
                       featureData = fd
                       )

#########################
#构建细胞发育轨迹
#########################
ordering_genes <- pbmc2[["RNA"]]@var.features
HSMM <- setOrderingFilter(HSMM, ordering_genes)
HSMM <- estimateSizeFactors(HSMM)
#降维
HSMM <- reduceDimension(HSMM,
                        norm_method="none", 
                        reduction_method="DDRTree",
                        max_components=3,
                        scaling=TRUE,
                        verbose=TRUE,
                        pseudo_expr=0)

#将细胞摆放到轨迹上
HSMM <- orderCells(HSMM)

#绘图展示
#以state来展示
pdf('trajectory_state.pdf')
plot_cell_trajectory(HSMM)
dev.off()

#以细胞类型来展示
pdf('trajectory_celltype.pdf')
plot_cell_trajectory(HSMM, color_by = 'celltype') 
dev.off()

#分别展示每一个细胞类型的发育轨迹
pdf('trajectory_celltype_separate.pdf')
plot_cell_trajectory(HSMM, color_by = 'celltype')+ facet_wrap(~celltype, nrow = 2) + NoLegend()
dev.off()

pdf('trajectory_celltype_separate2.pdf')
plot_cell_trajectory(HSMM, color_by = 'celltype')+ facet_wrap(~celltype, nrow = 2) 
dev.off()

#查看特定的gene在发育轨迹上的表达情况
pdf("5gene_trajectory.pdf",width = 20,height =20)
plot_cell_trajectory(HSMM,markers=c("INCENP","NEK11","AURKB","CCNA2","PRKAA2","CHEK1
","CDK5","PRKCD","CDK1"),use_color_gradient=T)
dev.off()

siggene<-c("INCENP","NEK11","AURKB","CCNA2","PRKAA2","CHEK1","CDK5","PRKCD","CDK1")

#heatmap
top5=scRNA.markers %>% group_by(cluster) %>% top_n(n = 5, wt = avg_log2FC)
sig_gene_names <- unique(top5$gene)
pseudotemporalplot<- plot_pseudotime_heatmap(HSMM[siggene,],
                                             num_clusters = 2,  #亚群数需要对应修改
                                             cores = 4,
                                             hmcols = NULL,
                                             show_rownames = T,
                                             return_heatmap = T)
pdf(file="monocle_heatmap.pdf")
pseudotemporalplot
dev.off()

pbmc2 <-subset(x =pbmc,celltype %in% c("Hepatocyte", 'Cancer stem cell'))

plot_cell_trajectory(HSMM, markers = c("VCAN","CD163","MRC1"),use_color_gradient = T,show_branch_points = F)

pseudotemporalplot<- plot_pseudotime_heatmap(HSMM[siggene,],
                                             num_clusters = 2,  #亚群数需要对应修改
                                             cores = 4,
                                             hmcols = colorRampPalette(rev(brewer.pal(9, "PRGn")))(62),
                                             show_rownames = T,
                                             return_heatmap = T)
pdf(file="monocle_heatmap2.pdf")
pseudotemporalplot
dev.off()


# 利用R包irGSEA (v2.1.5)进行单细胞功能富集分析，基于ssGSEA算法，计算GABA介导的巨噬细胞亚型细胞中HALLMARK通路得分，绘制得分热图进行展示